{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evolution\n",
    "\n",
    "Code examples from [Think Complexity, 2nd edition](http://greenteapress.com/wp/complexity2), Chapter 11\n",
    "\n",
    "Copyright 2016 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import itertools\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from thinkstats2 import Cdf\n",
    "from thinkstats2 import RandomSeed\n",
    "\n",
    "import thinkplot\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('animation', html='html5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Fitness landscape\n",
    "\n",
    "A genotype is represented by a length-N array of 0s and 1.\n",
    "\n",
    "The fitness landscape maps from each location in N-D space to a random fitness.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FitnessLandscape:\n",
    "    def __init__(self, N):\n",
    "        \"\"\"Create a fitness landscape.\n",
    "        \n",
    "        N: number of dimensions\n",
    "        \"\"\"\n",
    "        self.N = N\n",
    "        self.set_values()\n",
    "        \n",
    "    def set_values(self):\n",
    "        self.one_values = np.random.random(self.N)\n",
    "        self.zero_values = np.random.random(self.N)\n",
    "\n",
    "    def random_loc(self):\n",
    "        \"\"\"Choose a random location.\"\"\"\n",
    "        # in NumPy versions prior to 1.11, you might need\n",
    "        #return np.random.randint(2, size=self.N).astype(np.int8)\n",
    "        return np.random.randint(2, size=self.N, dtype=np.int8)\n",
    "    \n",
    "    def fitness(self, loc):\n",
    "        \"\"\"Evaluates the fitness of a location.\n",
    "        \n",
    "        loc: array of N 0s and 1s\n",
    "        \n",
    "        returns: float fitness\n",
    "        \"\"\"\n",
    "        fs = np.where(loc, self.one_values, self.zero_values)\n",
    "        return fs.mean()\n",
    "    \n",
    "    def distance(self, loc1, loc2):\n",
    "        return np.sum(np.logical_xor(loc1, loc2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, here's a 3-D landscape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fit_land = FitnessLandscape(3)\n",
    "fit_land.N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`one_values` and `zero_values` contain the fitness contributions of having a 1 or 0 at each element of the location array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_land.one_values, fit_land.zero_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitness of a location is the mean of its fitness contributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc = fit_land.random_loc()\n",
    "loc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.where(loc, fit_land.one_values, fit_land.zero_values)\n",
    "a, np.mean(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`fitness` evaluates the fitness of a location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc, fit_land.fitness(loc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`distance` computes the number of bit flips to get from one location to another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc1 = fit_land.random_loc()\n",
    "loc2 = fit_land.random_loc()\n",
    "print(loc1)\n",
    "print(loc2)\n",
    "fit_land.distance(loc1, loc2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It uses `np.logical_xor`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logical_xor(loc1, loc2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The agents\n",
    "\n",
    "Here's the class that represents agents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Agent:\n",
    "    \"\"\"Represents an agent in an NK model.\"\"\"\n",
    "    \n",
    "    def __init__(self, loc, fit_land):\n",
    "        \"\"\"Create an agent at the given location.\n",
    "        \n",
    "        loc: array of N 0s and 1s\n",
    "        fit_land: reference to an fit_land\n",
    "        \"\"\"\n",
    "        self.loc = loc\n",
    "        self.fit_land = fit_land\n",
    "        self.fitness = fit_land.fitness(self.loc)\n",
    "        \n",
    "    def copy(self):\n",
    "        return Agent(self.loc, self.fit_land)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each agent has a location, a reference to a FitnessLandscape, and a fitness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc = fit_land.random_loc()\n",
    "agent = Agent(loc, fit_land)\n",
    "agent.loc, agent.fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Simulator\n",
    "\n",
    "The `Simulator` class provides methods to run the simulations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Simulation:\n",
    "    \n",
    "    def __init__(self, fit_land, agents):\n",
    "        \"\"\"Create the simulation:\n",
    "        \n",
    "        fit_land: fit_land\n",
    "        num_agents: int number of agents\n",
    "        agent_maker: function that makes agents\n",
    "        \"\"\"\n",
    "        self.fit_land = fit_land\n",
    "        self.agents = np.asarray(agents)\n",
    "        self.instruments = []\n",
    "        \n",
    "    def add_instrument(self, instrument):\n",
    "        \"\"\"Adds an instrument to the list.\n",
    "        \n",
    "        instrument: Instrument object\n",
    "        \"\"\"\n",
    "        self.instruments.append(instrument)\n",
    "        \n",
    "    def plot(self, index, *args, **kwargs):\n",
    "        \"\"\"Plot the results from the indicated instrument.\n",
    "        \"\"\"\n",
    "        self.instruments[index].plot(*args, **kwargs)\n",
    "        \n",
    "    def run(self, num_steps=500):\n",
    "        \"\"\"Run the given number of steps.\n",
    "        \n",
    "        num_steps: integer\n",
    "        \"\"\"\n",
    "        # initialize any instruments before starting\n",
    "        self.update_instruments()\n",
    "        \n",
    "        for _ in range(num_steps):\n",
    "            self.step()\n",
    "        \n",
    "    def step(self):\n",
    "        \"\"\"Simulate a time step and update the instruments.\n",
    "        \"\"\"\n",
    "        n = len(self.agents)\n",
    "        fits = self.get_fitnesses()\n",
    "        \n",
    "        # see who dies\n",
    "        index_dead = self.choose_dead(fits)\n",
    "        num_dead = len(index_dead)\n",
    "        \n",
    "        # replace the dead with copies of the living\n",
    "        replacements = self.choose_replacements(num_dead, fits)\n",
    "        self.agents[index_dead] = replacements\n",
    "\n",
    "        # update any instruments\n",
    "        self.update_instruments()\n",
    "        \n",
    "    def update_instruments(self):\n",
    "        for instrument in self.instruments:\n",
    "            instrument.update(self)\n",
    "            \n",
    "    def get_locs(self):\n",
    "        \"\"\"Returns a list of agent locations.\"\"\"\n",
    "        return [tuple(agent.loc) for agent in self.agents]\n",
    "    \n",
    "    def get_fitnesses(self):\n",
    "        \"\"\"Returns an array of agent fitnesses.\"\"\"\n",
    "        fits = [agent.fitness for agent in self.agents]\n",
    "        return np.array(fits)\n",
    "    \n",
    "    def choose_dead(self, ps):\n",
    "        \"\"\"Choose which agents die in the next timestep.\n",
    "        \n",
    "        ps: probability of survival for each agent\n",
    "        \n",
    "        returns: indices of the chosen ones\n",
    "        \"\"\"\n",
    "        n = len(self.agents)\n",
    "        is_dead = np.random.random(n) < 0.1\n",
    "        index_dead = np.nonzero(is_dead)[0]\n",
    "        return index_dead\n",
    "        \n",
    "    def choose_replacements(self, n, weights):\n",
    "        \"\"\"Choose which agents reproduce in the next timestep.\n",
    "        \n",
    "        n: number of choices\n",
    "        weights: array of weights\n",
    "        \n",
    "        returns: sequence of Agent objects\n",
    "        \"\"\"\n",
    "        agents = np.random.choice(self.agents, size=n, replace=True)\n",
    "        replacements = [agent.copy() for agent in agents]\n",
    "        return replacements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use a few functions to create agents.  If we want to start with identical agents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_identical_agents(fit_land, num_agents, agent_maker):\n",
    "    \"\"\"Make an array of Agents.\n",
    "    \n",
    "    fit_land: FitnessLandscape\n",
    "    num_agents: integer\n",
    "    agent_maker: class used to make Agent\n",
    "    \n",
    "    returns: array of Agents\n",
    "    \"\"\"\n",
    "    loc = fit_land.random_loc()\n",
    "    agents = [agent_maker(loc, fit_land) for _ in range(num_agents)]\n",
    "    return np.array(agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or agents at random locations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_random_agents(fit_land, num_agents, agent_maker):\n",
    "    \"\"\"Make an array of Agents.\n",
    "    \n",
    "    fit_land: FitnessLandscape\n",
    "    num_agents: integer\n",
    "    agent_maker: class used to make Agent\n",
    "    \n",
    "    returns: array of Agents\n",
    "    \"\"\"\n",
    "    locs = [fit_land.random_loc() for _ in range(num_agents)]\n",
    "    agents = [agent_maker(loc, fit_land) for loc in locs]\n",
    "    return np.array(agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or one agent at each possible location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_all_agents(fit_land, agent_maker):\n",
    "    \"\"\"Make an array of Agents.\n",
    "    \n",
    "    fit_land: FitnessLandscape\n",
    "    agent_maker: class used to make Agent\n",
    "    \n",
    "    returns: array of Agents\n",
    "    \"\"\"\n",
    "    N = fit_land.N\n",
    "    locations = itertools.product([0, 1], repeat=N)\n",
    "    agents = [agent_maker(loc, fit_land) for loc in locations]\n",
    "    return np.array(agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`make_all_agents` uses `itertools.product`, which returns a generator that enumerates the Cartesian product of the set `{0, 1}` with itself `N` times, which is a fancy way to say that it enumerates all sequences of `N` bits.  Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_land = FitnessLandscape(3)\n",
    "agents = make_all_agents(fit_land, Agent)\n",
    "for agent in agents:\n",
    "    print(agent.loc, agent.fitness)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The distribution of fitness\n",
    "\n",
    "Let's create a fitness landscape with an agent at each location and see what the distribution of fitness looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "N = 8\n",
    "fit_land = FitnessLandscape(N)\n",
    "agents = make_all_agents(fit_land, Agent)\n",
    "sim = Simulation(fit_land, agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`plot_fitnesses` plots the CDF of fitness across the population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fitnesses(sim):\n",
    "    \"\"\"Plot the CDF of fitnesses.\n",
    "    \n",
    "    sim: Simulation object\n",
    "    \"\"\"\n",
    "    fits = sim.get_fitnesses()\n",
    "    cdf_fitness = Cdf(fits)\n",
    "    thinkplot.Cdf(cdf_fitness)\n",
    "    return np.mean(fits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initially the distribution is approximately Gaussian, because it's the sum of 8 independent uniformly distributed variates.  See the Central Limit Theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_fitnesses(sim)\n",
    "thinkplot.Config(xlabel='Fitness', ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After one time step, there's not much change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.step()\n",
    "plot_fitnesses(sim)\n",
    "thinkplot.Config(xlabel='Fitness', ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After 100 time steps, we can see that the number of unique values has decreased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.run(100)\n",
    "plot_fitnesses(sim)\n",
    "thinkplot.Config(xlabel='Fitness', ylabel='CDF')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instruments\n",
    "\n",
    "To measure these changes over the course of the simulations, we'll use Instrument objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Instrument:\n",
    "    \"\"\"Computes a metric at each timestep.\"\"\"\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.metrics = []\n",
    "        \n",
    "    def update(self, sim):\n",
    "        \"\"\"Compute the current metric.\n",
    "        \n",
    "        Appends to self.metrics.\n",
    "        \n",
    "        sim: Simulation object\n",
    "        \"\"\"\n",
    "        # child classes should implement this method\n",
    "        pass\n",
    "        \n",
    "    def plot(self, **options):\n",
    "        thinkplot.plot(self.metrics, **options)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MeanFitness` instrument computes the mean fitness after each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MeanFitness(Instrument):\n",
    "    \"\"\"Computes mean fitness at each timestep.\"\"\"\n",
    "    label = 'Mean Fitness'\n",
    "    \n",
    "    def update(self, sim):\n",
    "        mean = np.nanmean(sim.get_fitnesses())\n",
    "        self.metrics.append(mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's mean fitness as a function of (simulated) time for a single run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "N = 8\n",
    "fit_land = FitnessLandscape(N)\n",
    "agents = make_all_agents(fit_land, Agent)\n",
    "\n",
    "sim = Simulation(fit_land, agents)\n",
    "instrument = MeanFitness()\n",
    "sim.add_instrument(instrument)\n",
    "sim.run(500)\n",
    "sim.plot(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get a better sense of average behavior, and variation around the average, but plotting multiple runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sims(fit_land, agent_maker, sim_maker, instrument_maker, **plot_options):\n",
    "    \"\"\"Runs simulations and plots metrics.\n",
    "    \n",
    "    fit_land: FitnessLandscape\n",
    "    agent_maker: function that makes an array of Agents\n",
    "    sim_maker: function that makes a Simulation\n",
    "    instrument_maker: function that makes an instrument\n",
    "    plot_options: passed along to plot\n",
    "    \"\"\"\n",
    "    plot_options['alpha'] = 0.3\n",
    "\n",
    "    for _ in range(10):\n",
    "        agents = agent_maker(fit_land)\n",
    "        sim = sim_maker(fit_land, agents)\n",
    "        instrument = instrument_maker()\n",
    "        sim.add_instrument(instrument)\n",
    "        sim.run()\n",
    "        sim.plot(0, **plot_options)\n",
    "    thinkplot.Config(xlabel='Time', ylabel=instrument.label)\n",
    "    return sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`agent_maker1` puts one agent at each location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def agent_maker1(fit_land):\n",
    "    return make_all_agents(fit_land, Agent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With no differential survival or reproduction, we get a random walk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, Simulation, MeanFitness, color='blue')\n",
    "thinkplot.Save('chap11-1', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential survival\n",
    "\n",
    "We can add differential survival by overriding `choose_dead`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SimWithDiffSurvival(Simulation):\n",
    "    \n",
    "    def choose_dead(self, ps):\n",
    "        \"\"\"Choose which agents die in the next timestep.\n",
    "        \n",
    "        ps: probability of survival for each agent\n",
    "        \n",
    "        returns: indices of the chosen ones\n",
    "        \"\"\"\n",
    "        n = len(self.agents)\n",
    "        is_dead = np.random.random(n) > ps\n",
    "        index_dead = np.nonzero(is_dead)[0]\n",
    "        return index_dead"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With differential survival, mean fitness increases and then levels off."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithDiffSurvival, MeanFitness, color='blue')\n",
    "thinkplot.Save('chap11-2', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add differential reproduction by overriding `choose_replacements`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SimWithDiffReproduction(Simulation):\n",
    "\n",
    "    def choose_replacements(self, n, weights):\n",
    "        \"\"\"Choose which agents reproduce in the next timestep.\n",
    "        \n",
    "        n: number of choices\n",
    "        weights: array of weights\n",
    "        \n",
    "        returns: sequence of Agent objects\n",
    "        \"\"\"\n",
    "        p = weights / np.sum(weights)\n",
    "        agents = np.random.choice(self.agents, size=n, replace=True, p=p)\n",
    "        replacements = [agent.copy() for agent in agents]\n",
    "        return replacements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With differential reproduction (but not survival), mean fitness increases and then levels off."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithDiffReproduction, MeanFitness, color='blue')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** What if you have both?  Write a class called `SimWithBoth` that uses the new versions of `choose_dead` and `choose_replacements`.  Does mean fitness increase more quickly?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "class SimWithBoth(Simulation):\n",
    "    choose_dead = SimWithDiffSurvival.choose_dead\n",
    "    choose_replacements = SimWithDiffReproduction.choose_replacements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithBoth, MeanFitness, color='blue')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of different agents\n",
    "\n",
    "Without mutation, we have no way to add diversity.  The number of occupied locations goes down over time.\n",
    "\n",
    "`OccupiedLocations` is an instrument that counts the number of occupied locations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class OccupiedLocations(Instrument):\n",
    "    label = 'Occupied Locations'\n",
    "\n",
    "    def update(self, sim):\n",
    "        uniq_agents = len(set(sim.get_locs()))\n",
    "        self.metrics.append(uniq_agents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what that looks like with no differential survival or reproduction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, Simulation, OccupiedLocations, color='green')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** What effect do differential survival and reproduction have on the number of occupied locations?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithDiffSurvival, OccupiedLocations, color='green')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithDiffReproduction, OccupiedLocations, color='green')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "RandomSeed(17)\n",
    "plot_sims(fit_land, agent_maker1, SimWithBoth, OccupiedLocations, color='green')\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model we have so far might explain changes in existing populations, but it doesn't explain increasing diversity or complexity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mutation\n",
    "\n",
    "Mutation is one way of increasing, or at least maintaining, diversity.\n",
    "\n",
    "`Mutant` is a kind of agent that overrides `copy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Mutant(Agent):\n",
    "    \n",
    "    prob_mutate = 0.05\n",
    "    \n",
    "    def copy(self):\n",
    "        if np.random.random() > self.prob_mutate:\n",
    "            loc = self.loc.copy()\n",
    "        else:\n",
    "            direction = np.random.randint(self.fit_land.N)\n",
    "            loc = self.mutate(direction)\n",
    "        return Mutant(loc, self.fit_land)\n",
    "    \n",
    "    def mutate(self, direction):\n",
    "        \"\"\"Computes the location in the given direction.\n",
    "        \n",
    "        Result differs from the current location along the given axis.\n",
    "        \n",
    "        direction: int index from 0 to N-1\n",
    "        \n",
    "        returns: new array of N 0s and 1s\n",
    "        \"\"\"\n",
    "        new_loc = self.loc.copy()\n",
    "        new_loc[direction] ^= 1\n",
    "        return new_loc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test it out, I'll create an agent at a random location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 8\n",
    "fit_land = FitnessLandscape(N)\n",
    "loc = fit_land.random_loc()\n",
    "agent = Mutant(loc, fit_land)\n",
    "agent.loc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we make 20 copies, we expect about one mutant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(20):\n",
    "    copy = agent.copy()\n",
    "    print(fit_land.distance(agent.loc, copy.loc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`agent_maker2` makes identical agents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def agent_maker2(fit_land):\n",
    "    agents = make_identical_agents(fit_land, 100, Mutant)\n",
    "    return agents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we start with identical mutants, we still see increasing fitness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "sim = plot_sims(fit_land, agent_maker2, SimWithBoth, MeanFitness, color='blue')\n",
    "thinkplot.Save('chap11-3', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now the number of occupied locations increases, reaching a steady state at about 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "sim = plot_sims(fit_land, agent_maker2, SimWithBoth, OccupiedLocations, color='green')\n",
    "thinkplot.Save('chap11-4', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In steady state, many agents are at the optimal location, and others are usually just a few mutations away.  To quantify that, we can compute the mean distance between all pairs of agents.\n",
    "\n",
    "The distance between two agents is the number of bit flips to get from one location to another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MeanDistance(Instrument):\n",
    "    \"\"\"Computes mean distance between pairs at each timestep.\"\"\"\n",
    "    label = 'Mean Distance'\n",
    "        \n",
    "    def update(self, sim):\n",
    "        N = sim.fit_land.N\n",
    "        i1, i2 = np.triu_indices(N)\n",
    "        agents = zip(sim.agents[i1], sim.agents[i2])\n",
    "        \n",
    "        distances = [fit_land.distance(a1.loc, a2.loc)\n",
    "                     for a1, a2 in agents if a1 != a2]\n",
    "        \n",
    "        mean = np.mean(distances)\n",
    "        self.metrics.append(mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean distance is initially 0, when all agents are identical.  It increases as the population migrates toward the optimal location, then settles into a steady state around 1.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "fit_land = FitnessLandscape(10)\n",
    "agents = make_identical_agents(fit_land, 100, Mutant)\n",
    "sim = SimWithBoth(fit_land, agents)\n",
    "sim.add_instrument(MeanDistance())\n",
    "sim.run(500)\n",
    "sim.plot(0, color='orange')\n",
    "thinkplot.Config(xlabel='Time', ylabel='Mean distance')\n",
    "thinkplot.Save('chap11-5', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changing landscape\n",
    "\n",
    "One cause of speciation is a change in the landscape.  Suppose a population in steady state in one landscape is transported to another landscape.  After a period of migration, it would settle around the new equilibrium point, with a small mean distance between agents.  The mean distance between the original agents and the migrants is generally much larger.\n",
    "\n",
    "The following simulation runs 500 steps on one landscape, then switches to a different landscape and resumes the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RandomSeed(17)\n",
    "\n",
    "fit_land = FitnessLandscape(10)\n",
    "agents = make_identical_agents(fit_land, 100, Mutant)\n",
    "sim = SimWithBoth(fit_land, agents)\n",
    "sim.add_instrument(MeanFitness())\n",
    "sim.add_instrument(OccupiedLocations())\n",
    "sim.add_instrument(MeanDistance())\n",
    "sim.run(500)\n",
    "locs_before = sim.get_locs()\n",
    "fit_land.set_values()\n",
    "sim.run(500)\n",
    "locs_after = sim.get_locs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the switch to a new landscape, mean fitness drops briefly, then the population migrates to the new optimal location, which happens to be higher, in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot(0, color='blue')\n",
    "plt.axvline(500, color='gray', linewidth=3, alpha=0.5)\n",
    "thinkplot.Config(xlabel='Time', ylabel='Mean Fitness')\n",
    "thinkplot.Save('chap11-6', clf=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of occupied locations (sometimes) increases while the population is migrating."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot(1, color='green')\n",
    "plt.axvline(500, color='gray', linewidth=3, alpha=0.5)\n",
    "thinkplot.Config(xlabel='Time', ylabel='Occupied Locations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the mean distance (sometimes) increases until the population reaches the new steady state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.plot(2, color='orange')\n",
    "plt.axvline(500, color='gray', linewidth=3, alpha=0.5)\n",
    "thinkplot.Config(xlabel='Time', ylabel='Mean Distance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean distance between clusters is much bigger than the dispersion within clusters, so we can interpret the clusters as distinct species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = []\n",
    "for loc1 in locs_before:\n",
    "    for loc2 in locs_after:\n",
    "        distances.append(fit_land.distance(loc1, loc2))\n",
    "np.mean(distances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** When we change the landscape, the number of occupied locations and the mean distance usually increase, but the effect is not always big enough to be obvious.  You might want to try out some different random seeds to see how general the effect is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
